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1. Introduction 

CP violation is now a well established property of Nature, having been observed in both kaon 
and B meson systems. Its standard model explanation, that the violation arises from the difference 
between the mass and weak interaction eigenstates of the quarks, is quantitatively encapsulated 
in the unitary 3x3 Cabbibo-Kobayashi-Maskawa (CKM) matrix which relates these eigenstates. 
Determining the elements of this matrix is the subject of extensive theoretical and experimental 
effort. However, even assuming precise values for the elements of the CKM matrix, a prediction 
for the rate of direct CP-violating decays in kaons, determined experimentally by measuring the 
quantity Re(e'/e), is lacking. Such a prediction relies on values for the hadronic K — > nn matrix 
elements of low-energy, four-quark operators. Lattice QCD is the only first principles method to 
allow determination of these matrix elements and their measurement has been the subject of many 
studies. In this note, we discuss our progress in determining these matrix elements, using full 2+1 
flavor QCD with domain wall quarks. 

Euclidean space lattice methods for the direct evaluation of the desired K — > nn matrix ele- 
ments have been developed, and some preliminary studies have been performed [1,2]. However, to 
reach physical values for the kaon and pion masses in full QCD simulations is beyond the reach of 
current computers. Here we follow the long-standing approach of using chiral perturbation theory 
to relate %% matrix elements to matrix elements of the same operator measured in K — > % and 
K — > matrix elements [3]. Quenched calculations have been done with this approach [4, 5], which 
showed that with DWF, the problems of lower dimensional operator mixing at finite lattice spacing 
and renormalization of the operators were under good control. Statistically well resolved values 
for the lowest order chiral perturbation theory constants were achieved and, by naively extending 
lowest order chiral perturbation theory to the kaon mass, physical results were quoted. 

However, for these particular kaon weak matrix elements calculations, it was pointed out that 
quenching makes the relation between the full QCD chiral perturbation theory constants and the 
quenched ones ambiguous [6, 7, 8], providing a source of uncertainty in interpreting the quenched 
calculation that is more serious than the general caution that must be exercised with quenching. In 
the current 2+1 flavor calculation, these ambiguities are removed. 

To date, our calculation involves measuring the same operators as in our previous quenched 
calculation [4], although on larger volume lattices, with 2+1 dynamical flavors and with much 
lighter valence quark masses [9, 10]. With these valence masses, we see that next-to-leading order 
(NLO) SU(3) chiral perturbation theory for pseudoscalar decay constants and masses [11], and 
the pseudoscalar bag parameter Bps (the analogue of Bk with two light quarks in the "kaon") 
[12], represents the lattice data at the « 5% level or better. Thus we can try similar NLO chiral 
perturbation theory fits to the data presented here, which covers a similar range of pseudoscalar 
masses, in the hopes of extracting the desired low energy constants. These fits require the 2+1 
flavor partially quenched chiral perturbation theory formula for our matrix elements, a calculation 
which is underway [13]. 

We also discuss the status of our calculation of the non-perturbative renormalization Z factors 
needed to convert our lattice results to MS conventions. 
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Figure 1: The contractions used in the evaluation of the AS = 1 operators that have been calculated in this 
report. The sources for the kaons and pions are Coulomb gauge fixed wall sources and the closed fermion 
loop contractions are calculated with stochastic sources. 

2. Simulation Details 

We have measured the 10 AS = 1 operators needed for a determination of e' / £ and the AI = 1/2 
rule, using the notation of [4]. When evaluated between K + and 7i + states, the operators produce 
the contractions shown in Figure 1. The contractions between kaons and pions are known as figure 
eight diagrams and eye diagrams. We also need vacuum diagrams and contractions involving sd, to 
remove the mixing with lower dimensional operators. We use Coulomb gauge fixed wall sources at 
t = 5 and 59 on our 24 3 x 64 x 16 lattices as the sources for the kaons and pions. The Dirac equation 
is solved twice for each wall source, once with periodic and once with anti-periodic boundary 
conditions, and the two solutions added to effectively double the length of the lattice in the time 
direction. Our operator contractions are always done with 5 < t op < 59, where t op is the time where 
the operator is evaluated, so the double lattice merely serves to suppress any contributions to our 
operators from propagation around the world in the time direction. 

For the closed fermion loops in the eye and K — > contractions, we use a stochastic source, 
spread over the spatial volume and covering all time slices from 12 to 51 (40 time slices). A single 
stochastic source is a fixed spin and color component, with a space-time value that is Gaussianly 
distributed. For each spin and color combination, solutions are found for four different stochastic 
space-time distributions. 

For each configuration, 6 valence quark masses are used: nif = 0.001, 0.005, 0.01, 0.02, 0.03, 
0.04. The residual quark mass from the finite L s extent of our lattices is 0.00315(2) [11]. These bare 
quark masses correspond to pseudoscalar masses ranging from 250 MeV to 750 MeV. From fits to 
pseudoscalar masses and decay constants, we find that for the three lightest bare quark masses 
(0.001, 0.005 and 0.01), NLO chiral perturbation theory agrees well with our measured data. 

For each configuration, we do 8 full (all spins and colors for the source) five-dimensional 
solutions to the Dirac equations for each valence mass. Each wall source requires 2 solutions (one 
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Irrep 


Number 


Isospin 


K + -> 7T+ 




(27,1) 


1 


1/2, 3/2 


_4< a (27,l) 


4/ 2 ^(27,1) 


(8,1) 


4 


1/2 




^m 2 a (8,1) 


(8,8) 


2 


1/2, 3/2 







Table 1: The relation between matrix elements and low energy constants, to lowest order in chiral perturba- 
tion theory. The second column gives the number of different operators transforming as the given irreducible 
representation. 



for each boundary condition) and we do four stochastic estimators for the closed fermion loops. 
Thus we spend half of the time on the wall source inversions and half on the stochastic inversions. 

The results reported here are for 76 configurations, separated by 40 molecular dynamics time 
units, for our ensemble with light dynamical quark mass m/ = 0.005 and 74 configurations, with 
the same separation, for m/ = 0.01. To do these 8 full inversions for the six valence quark masses 
on a single configuration takes about 24 hours on a 4,096 node QCDOC partition, with the lightest 
valence quark mass taking about half of the total time. 



3. Chiral Perturbation Theory and Ward Identity Test 

In the 3 flavor effective theory that comes from integrating out the charm, bottom and top 
quarks (as well as the weak vector bosons), the effective Hamiltonian for AS = 1 transitions is 

jftAS=i) = | u ; j|< [zf{jl) + Ty . m a J (3.!) 

Here z,-(jll) and y,-(ju) are Wilson coefficients and the Qi are four-quark operators. Of the ten four- 
quark operators <2, above, only seven are independent, and they include three representations of 
SU(3)l x SU(3)r: (27,1), (8,1) and (8,8). In lowest order chiral perturbation theory, the desired 
K — ► KK matrix elements are all proportional to an appropriate low energy constant, denoted as 

a[ S ' l) , at 27 - 1 ), or a( 8 - 8 ). These come from K — ► n and K — ► matrix elements as given in Table 

(8 1) 

1. Note that for the (8,1) operators, to determine the physically relevant quantity a\ ' requires 

(8 1) 

remove the power divergent quantity ' , which is measured in K — > matrix elements 

To accurately determine the leading order chiral perturbation theory constants, we want to try 

to fit our matrix element data to the full partially quenched 2+1 flavor formulae. Calculations of 

these formulae are underway [13], and when they are complete, we will fit our data to them in our 

efforts to extract the best values for the a's. 

Since NLO chiral perturbation theory plays such an important role here, we turn briefly to a 

simple matrix element (n + \sd\K + ) that is known in chiral perturbation theory and that also enters 

our calculation. At leading order in chiral perturbation theory, one has 

9/77 4-' 

—f(n + \sd\K + )-\ =0 (3.2) 
mi 
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Q-O N f = 0, a" = 1.92 GeV 

□— □ N f = 2+1, m, = 0.005, a"' = 1.73 GeV 

0—0 N f = 2+1, m, = 0.01, a"' = 1.73 GeV 
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Figure 2: A plot of the appropriately scaled (n + \sd\K + ) versus valence quark mass. When the valence and 
sea quark masses are equal, the curves should go to 1 . 



We calculate this from the ratio of a three point function to two, 2-point functions as 

(P%f(xo)[sd](y)P^(zo)) 



R 7 



(P^(xo)P^{zo)) 



(3.3) 



Here P™ 11 is a wall Coulomb gauge fixed pseudoscalar wall source with the quantum numbers of 
the n + . /?2 differs from (n + \sd\K + ) by known factors. 

Figure 2 shows a plot of the appropriately scaled R2 versus valence input quark mass for the 
current 2+1 flavor full QCD ensembles and the same results for our earlier quenched ensemble. In 
the limit m/ — ► — m res and the light sea and valence quarks are equal, the curves should go to 1. 
We can see that there is noticeable curvature in the graphs for quark masses below 0.01 and there 
is a small dynamical quark mass effect when comparing the m/ = 0.005 and 0.01 ensembles. This 
is certainly the region where partially quenched chiral logarithms should be noticeable, and are 
seen in our other observables, but fits to the precise formulae (when available) will be required to 
sharpen this general statement. For Bps [12] we do see that the partially quenched logarithms have 
a curvature that is opposite the curvature when the valence and sea quark masses are the same. Also 
note that the curvature is much more noticeable in the dynamical ensemble data than the quenched 
results and that the valence masses in the current work are much lighter than the earlier quenched 
work. 



4. AI = 3/2 Matrix Elements 

We now turn to the AI = 3/2 matrix elements and present some of our lattice data. The 3/2 
part of Q2 and Q& are presented in Figure 3. We see at once that we have very long plateaus which 
extend, with essentially uniform errors, across the lattice. We choose to fit the plateau from t = 12 
to 51, i.e. starting from a distance of 7 away from the two wall sources. Clearly with such long 



5 



Kaon Weak Matrix Elements in 2+1 flavor DWF QCD 



Robert D. Mawhinney 



0.0015 



\4 



o 



0.001 - 



0.0005 




-0.02 



-0.03 



-0.04 



-0.05 



-0.06 




(a) 



10 20 30 40 50 
t 

(b) 



Figure 3: Plateaus for the AI = 3/2 part of the K — ► % matrix elements of Q2 and Q%. For Q2 (left panel) the 
valence masses run from the largest (top points) to the smallest (bottom), while for Q% (right panel) they run 
from the smallest (top points) to the largest (bottom). The pseudoscalars have degenerate quark masses. 



plateaus the results are not sensitive to the precise fit range chosen. The average value, with error 

bars, for each mass is indicated by the horizontal lines on the graphs. 

(3 2) 

In Figure 4, we present values for Q 2 as a function of valence quark mass for our previous 
quenched calculation (left panel) and the mi = 0.005 lattices from our current ensemble. Note the 
much smaller range of quark masses plotted in the 2+ 1 flavor QCD case. In our earlier quenched 
work, we saw no strong indication of quenched chiral logarithms, but found a large contribution 
from the analytic NLO terms. We do not yet have the chiral perturbation theory formulae for the 
2+1 flavor case, nor the fits to the data. However, one can see that some downward curvature is 
needed to have the extrapolated matrix element vanish at the origin, as required by chiral symmetry. 
Also, for these much smaller valence quark masses, one does not readily see a noticeable effect of 
higher order analytic terms. 



5. Weak Matrix Elements in the (8,1) Representation 
5.1 (8,1) Operators 

Among the 10 weak operators involved in the K — > %% process below the charm threshold, 
Qi {/ = 3,4,5,6} are pure (8,1) AI = 1/2 operators, while Q t {i = 1,2,9,10} have both (8,1), AI = 
1 /2 and (27,1), AI = 1/2 parts. The (8,1) operators are particularly important for the investigation 
of CP violation, since, from the size of the Wilson coefficients at 2 GeV, the matrix element of 
is dominant in Im (Aq) [4]. 

Since the AI = 1/2 operators can mix with the quadratically divergent operator @( 3 ^) = 
s (1 — 75) d, their K — > 71 matrix elements alone cannot determine the relevant K — > %% elements. It 
is necessary to evaluate their K —> matrix elements to remove the unwanted contributions. 
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Figure 4: A comparison of the lattice normalized matrix elements for our earlier quenched calculation (left 
panel) and the current 2+1 flavor calculation (right panel). 



5.2 Plateaus of K — » % for (8,1) Operators 

As mentioned earlier, we evaluate the value of (n + \Qj\ K + ), by placing Coulomb gauge fixed 
wall sources at tx = 5 and t % = 59. The eye contractions that appear in (8,1) operators are deter- 
mined with random sources with support from t ranc i = 12 to 51. We can then evaluate the operator 
contribution from any timeslice given by t op with 12 < t op < 51. By taking the ratios of the three- 
and two-point pseudoscalar-pseudoscalar Green's functions we calculate the required matrix ele- 
ments [4]. As an example, Figure 5 shows the plateaus for (n + |<2e|^ + )- With either degenerate 
or non-degenerate valence masses, we have very long plateaus (40 time-slices). Comparing these 
graphs with the plateaus in Figure 3, we see that the error bars here, after averaging over time slices 
(black horizontal lines) are smaller than the error bars for the individual points. Such an effect is 
not very noticeable in the AI = 3/2 amplitudes of Figure 3. 

5.3 Resolving Quadratic Divergence with K —> Matrix Elements 

As explained earlier, we need to calculate the K — > matrix elements for AI = 1/2 operators 
to remove the mixing with lower dimensional operators. For the (8,1) operators, we determine a 
mixing coefficient Tji ( as in [4] by calculating 





n (l/2) 
t£i',lat 




(0|( 





r\Q,i + r\\,i{m s -m d ) (5.1) 

where the valence quark masses have m s ^ ma- Again using Q§ as an example, Figure 6 shows the 
plateau of the K — > matrix element in the left panel and the result of fitting to Eq. 5. 1 above. One 
sees that the data shows a very accurate linear dependence on m s — nid, which is expected since the 
quadratically divergent part of dominates the numerator and it is proportional to m s — m c {. Since 
<26 has the largest divergent contribution, its matrix elements show the best linearity. For other 
operators, with a smaller divergent part, the linearity is less good, but also the required cancellation 
in the subtraction of the divergent part is less stringent. 
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Figure 5: Plateaus of (n + \Q(,\K + ) with m sea = 0.005. Figure (a) has degenerate valence masses (mj = 
m s = 0.001) and Figure b) has non-degenerate valence quark {m c [ = 0.001, m s = 0.005). In both cases, the 
matrix elements have long plateaus. 
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Figure 6: Calculating mixing coefficient between four-quark operators and sd using the ratio 



0(1/2) 



K°J I (0|(575c/) lat |^r°). Figure (a) shows the plateau of the matrix element with m sea = 0.005, 
= 0.001, and m s = 0.005. Figure (b) shows the results of fitting to (m s — m c [), which reveals a nice linear 
structure. 



5.4 Subtraction of the Quadratic Divergence 



After calculating the mixing coefficient T]i ,-, we use it to subtract the quadratically divergent 
contribution of (n + \sd\ K + ) from the raw matrix elements [4]. 



K 



or 



sub 



K 



4a 



Qf }) K + ) + T] u {m s + m d ) (n + \sd\K* 



f 2 



-mKtn n + (chiral logs) + (analytic terms) 



(5.2) 
(5.3) 
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Figure 7: Subtraction of quadratic divergence from/f — > % matrix elements. Figure (a) shows the subtraction 
process for degenerate valence masses mj = m s . It is clear that for the divergent operator, the physical 
quantity is very sensitive to any noise in the matrix elements. Figure (b) fits the subtracted matrix elements 
Ifi6 1 K + ) sub with all available pairs of non-degenerate valence masses to the leading order term mK>n K . 
The upper line does not use m res while on the lower line each subtraction has used m res . Both figures present 
the data of m sea = 0.005 ensemble. 

where the second line is the prediction of the chiral perturbation theory. (The cancellation of 
the power divergences is true without reference to chiral perturbation theory.) While we have 
begun fitting to next-to-leading order 2+ 1 flavor partially quenched chiral perturbation theory, the 
formulae are still being checked, so we have also experimented with simple leading order fits. 

Figure 7(a) shows the subtraction process for <26- Since is the most divergent, the sub- 
traction term (green squares) is very close to the unsubtracted matrix elements, and thus a precise 
determination of the subtraction term is vital. (This was also the case in the quenched calcu- 
lations.) In Figure 7(b), we show the subtracted matrix elements (ti + \Q$\ K + } sub , including all 
non-degenerate combinations of valence masses (mj, m s ), and also the result of fitting them to the 
leading order term in chiral perturbation theory, m^m^. The plots show that the data fit the lead- 
ing order theory very well. Of course, here our error bars are larger than for the AI = 3/2 matrix 
elements and modest chiral logarithm effects are consistent with the apparent linearity of our data 
and our error bars. It remains an open question as to whether we can achieve statistical accuracy 
for the subtracted (8,1) operators that is the same size as the predicted effects of chiral logarithms. 

Due to the existence of power divergences and the fact that the residual chiral symmetry break- 
ing effects in divergent amplitudes are not strictly proportional to m res , the subtracted matrix ele- 
ments do not have to vanish at the chiral limit [4]. (Recall that the physical quantity of interest is 
the slope of the subtracted matrix elements.) Rather, at the chiral limit, their deviation from zero 
should be an (? (m res ) effect. In Figure 7(b), the top line does not use m KS in the subtraction process 
(Eq.5.2), while the bottom line does subtraction with m sea . As these two lines bracket the zero point 
at the chiral limit, the above expectation is verified. 

5.5 PQS vs. PQN 

As we work in the 2+1 flavor partially quenched theory, there are different methods to make 
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Figure 8: Quark flow diagrams needed to evaluate the weak matrix elements for Qj {i = 3,4,5,6}. A 
connected line represents a trace over spin indices. Only diagram (b) has differs between the PQS and PQN 
method. With the PQS method, the quark loop (red) will be constructed from sea quarks, while with the 
PQN method, it will be constructed from valence quarks. 



Slope of <7r+|e,-|^ + ) sub 


3 


4 


5 


6 


PQS 


0.0009(37) 


0.0383(27) 


-0.0626(53) 


-0.183(13) 


PQN 


-0.0042(38) 


0.0373(34) 


-0.0638(55) 


-0.165(14) 



Table 2: Comparison between the PQS and PQN methods of the slope of the subtracted K — > n matrix 
elements (w.r.t m^m^). 



contractions with the four-quark operators which have a singlet piece, i.e. PQS vs. PQN [8, 14, 15, 
16]. The operators Qi {i = 3,4,5,6} all have a singlet part in their operator definition, 

Qi = {sd) L £ (qq) La[R . (5-4) 

qe{u,d,s} 

When evaluating their weak matrix elements on our partially quenched lattices, we may either keep 
the singlet structure of the operator, such that the summation will become a sum over all valence, 
sea, and ghost quarks, or we only sum over the valence quarks as we shall do in full QCD. The 
former, as we keep the singlet property, is called the PQS (partially quenched singlet) method, 
while the latter is called the PQN (partially quenched non-singlet) method. 

On the quark level, calculating the weak matrix elements of the above operators requires the 
evaluation of the three quark-flow diagrams shown in Figure 8. The only difference between PQS 
and PQN is in diagram (b). With the PQS method, the summation on the quark self-contraction (the 
red quark loop) will go over all valence, sea, and ghost quarks (and since valence loops cancel with 
ghosts, only the sum over sea remains). With the PQN method, the sum will go over the valence 
quarks only. 

We have calculated the subtracted matrix elements with both the PQS and PQN methods (Table 
2). The comparison of the slope of the subtracted matrix elements (w.r.t m^m^), which is the 
leading order coefficient for (8,1) operators, shows small differences between the two methods. 
We see that this ambiguity does not affect the chiral fits significantly on our 2+1 flavor partially 
quenched lattices. 

6. Non-perturbative Renormalization for AS = 1 Operators 
6.1 NPR Procedure 

To convert the bare lattice quantities into the continuum scheme, we have employed the Rome- 
Southampton non-perturbative renormalization prescription [17]. Following the procedure as de- 
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scribed in ref. [4], we write the NPR formula as 



or t, re n (Ai)= £ z .. (Ai) 



lat 



+ 0{a) 



(6.1) 



where 0, are the dimension 6 operators, and are the lower dimensional operators. The NPR cal- 
culation is done on 16 3 x 32 x 16 lattices [18] which have been generated with the same parameters 
as our 24 3 x 64 x 16 lattices, except that the lightest dynamical quark mass on the smaller volume 
is 0.01. 

To simplify the computation, we have eliminated Q4, Qg, and Q\o from the 10-operator basis, 
since they can be written as linear combination of the other operators [4]. Further, we rotate the 
remaining 7-operator basis such that each operator is in a distinct SU (3) (g> SU (3) representation. 
The transformation relation can be found in ref. [4]. The new set of operators is denoted as Q' t 

6.2 Resolving the Mixing with Lower Dimensional Operators 

To resolve the mixing between the four quark operator and the lower dimensional operators, 
we consider the two most divergent lower dimensional ones [4] : 

B\ = sd 

fi 2 = s(-0 +m^d + s{jj> + m^jd. (6.2) 
To compute the corresponding mixing coefficients c\ and c' 2 , we have used two conditions 



Tr 



Tr [(s(p)Or b d-(p)) amp 

if{s{ P )0T h d{ P )) amv 



= 
= 



(6.3) 



with propagators evaluated at unitary masses (m^ = m s = m sea ). 



6.3 Resolving the Mixing between Four-Quark Operators 

To determine the mixing coefficient between the 7 operators in the new basis, we construct a 
set of external quark combinations s [4]. Then we impose the renormalization condition that 

1 



Z 2 A 1 PaSy 



QSUb E J 



ajiyS 



amp 



(6.4) 



where ^ J p a5r is the projector corresponding to the spin and color structure of the operator j, and 
there is no sum over j in the above equation. On the r.h.s, F k i is the free field limit of the matrix 
^PaSy^f^^aPyd) • ^ ^ ^ s tne 1 uar k wavefunction renormalization factor [19]. 

The calculation of the combination r^ Q , 5)/ ^0- ub £'^^^ requires evaluation of the many 

quark flow diagrams that can possibly be constructed with the relevant operator and the external 
quark combination. Figure 9 shows 4 out of 248 quark flow diagrams involved in the evaluation 
of renormalization coefficient for Q' v Each diagram is evaluated at unitary quark mass and non- 
exceptional external momenta p\ and pi where p\ = p\ = (p\ — P2) [4]. 
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Figure 9: 4 out of 248 quark flow diagrams with their corresponding factors involved in the calculation 
of the mixing coefficients between Q\ and Q[ . The solid lines represent quark propagators and each loop 
formed by solid and dashed lines represents a trace over spin and color indices. 
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Figure 10: Evaluation of the mixing coefficients between four quark operators. Figure (a) shows the linear 



mass fit of the mixing coefficient Z q Z\ \. Figure (b) shows the final value of Z q Z\\ at 5 different momen- 
tum scales where we have sufficient pairs of non-exceptional momenta (pi, pj) to evaluate the quark flow 
diagrams. Also in Figure (b), we show a comparison between linear and quadratic mass fits, as well as a 
comparison between the values before and after the subtraction of c\. It is clear that the differences are not 
significant. 



After the evaluation of the entire matrix of mixing coefficients, we remove the coefficients 
that are both theoretically suppressed and statistically zero, and then fit the remaining matrix to a 
unitary mass dependence (Figure 10(a)). Since one of the lower dimensional mixing coefficients, 
c\, is expected to have the mass dependence m/a 2 at leading order, this step will automatically 
remove c\. Then the subtraction of the other, c\, is performed. Due to its small magnitude, this 
subtraction does not make big differences in the final result (Figure 10(b)). 



7. Conclusions 

We have discussed the status of measurements of K — > % and K — > weak matrix elements on 
24 x 64, 2+1 flavor dynamical domain wall lattices by the RBC and UKQCD collaborations. We 
are working at much smaller masses and larger volumes than earlier quenched calculations and, 
with the 2+1 flavor ensembles, we are free of quenching errors. We find good statistical errors and 
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long plateaus for our measured matrix elements and the subtraction of power divergent operator 
mixing contributions to (8,1) operators works as expected. We find that the ambiguity between 
the PQS and PQN methods is removed. The non-perturbative renormalization of our operators, 
following the Rome-Southampton prescription, is essentially complete. We await the completion 
of the NLO PQ%PT formulae for these matrix elements, and preliminary fits are underway. The 
accuracy and stability of these fits will play a large role in determining our final errors. 
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